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Abstract 

A new method for numerical solving of boundary problem for ordinary 
differential equations with slowly varying coefficients which is aimed at better 
representation of solutions in the regions of their rapid oscillations or exponential 
increasing (decreasing) is proposed. It is based on approximation of the solution to 
iind in the form of superposition of certain polynomial-exponential basic functions. The 
method is studied for the Heknholtz equation in comparison with the standard finite 
difference method. The numerical tests have shown the convergence of the method 
proposed. In comparison with the finite difference method the same accuracy is 
obtained on substantially rarer mesh. This advantage becomes more pronounced, if the 
solution varies very rapidly. 

1. Introduction 

In many fields of physics the problem of wave propagation and absorption in 
non-uniform media is under investigation. For example, in the controlled fusion 
problem the wave propagation phenomena are analyzed both analytically and 
numerically in application to the radio frequency plasma heating (e.g. [1]), the MHD 
plasma stability (e.g. [2]), the radio frequency plasma production (e.g. [3]). The 
problem of wave propagation is described, as a rule, by a set of differential equations. 
These equations are usually solved using the Fourier expansion or the discretization 



(e.g. [4]). In the case of rapidly oscillating or exponentially increasing (decreasing) 
wave fields, employing of the standard discretization methods requires a very fine 
mesh. For example, to simulate the radio frequency field excitation in plasma of the 
LHD steUarator [3] the mesh with number of nodes A^=20000 was used. This number 
will increase, if such modelling is performed in a larger device. 

There exist a lot of problems when the coefficients of the set of differential 
equations which describe the wave propagation problem vary much slower than the 
solutions. In this case, in the regions of rapid oscillations of the wave field one can use 
the WKB solutions (e.g. [5]). However, using them, one cannot provide prescribed 
accuracy and, therefore, cannot control convergence of the calculations. For this 
reason the WKB approximation cannot be considered as a proper numerical method. 

In the present paper we propose a new method for solving one-dimensional 
problem which improves the accuracy of numerical solution especially in the regions of 
rapid oscillations of the wave field. We call this approach as the local solution method. 
This method essentially exploits the fact that the solution of a system of Unear 
differential equations can be represented as superposition of some linearly independent 
solutions of the homogeneous system and a specific solution. The main idea is to 
approximate this representation using certain basic functions in every mesh cell. There 
is wide freedom in choice of such functions. When using them in the form of 
polynomials we obtain results similar to the results of the standard mesh methods. 
Aiming at more exact representation of the solution in the regions of its rapid 
oscillation and exponential increase (decrease) we have chosen the basic functions in a 
poljTiomial-exponential form suggested by the form of the WKB solutions. 

The method proposed is analyzed for the Helmholtz equation, the simplest 
equation describing wave processes. In section 2 we analyze theoretically three 
versions of the local solution method. In section 3 we present the results of numerical 
tests of the method in comparison with results obtained by the standard finite 
difference method. 

2. Formulation of the local solution method 

Consider the one-dimensional Helmholtz equation 



72 

-^y{x) + G{x)y{x) = R{x) , (1) 

which is defined at the interval x e (jc^jcJ. We assume that the function G{x) is not 

rapidly oscillating and has no breaks or singular points at this interval. To solve 
equation (1) numerically we introduce a mesh with the nodes x/ where i=l,2, ... ,n. 
The strategy of the numerical solving of this equation is as follows. Within the segment 
Si ( X, < X < x,^j ), the solution of equation (1) can be written in the form 

/■\x) = CfVr>(x) + C<'>y«(x) + /^\x) , (2) 
where yl'^{x) and y2^{x) are two linearly independent solutions of the homogeneous 
equation (R(x)-0) and y^^ix) is a specific solution. Making use of the smaUness of the 
segment considered, we assume that we can find approximations yf^x) (here j=l, 2, 
R) for these solutions with prescribed accuracy: 

yf{x)=yf{x) + 0(}r^), (3) 
where h=Xi+i-Xi, m is the degree of the approximation. Since the approximate solutions 
y]'\x) are known, in order to obtain the approximate solutions of equation (1), i.e. 

to determine the unknown coefficients Cf^ and 

rKx) = c^^'y^ix) + c^'y^Kx) + nKx) , (4) 

at every segment we have to match the solutions and their derivatives at the internal 
mesh nodes: 

dy^'^'' 
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dx 



dx 



(5) 



In this way we obtain 2(n-2) linear algebraic equations for 2(n-l) unknowns. Two 
equations stiU needed are to be obtained from two boundary conditions at the end 
point of the interval xi , Xr . The resulting matrbc of the described equation set has 
narrow band and, therefore, can be easily reversed. For the case of three segments and 
the boundary conditions posed at the opposite ends of the interval, the portrait of 
matrix is shown in equation 1. 
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Figure 1. The portrait of the matrix of the system. 

The simplest way to obtain the approximate solutions yj'\x) is to represent 

them in the form of power series. The unknown coefficients of the series could be 
obtained substituting these approximate solutions into equation (1) which is taken 
homogeneous to find functions y^'^ix) and y2^{x)and non-homogeneous for the 
specific solution y^^x), and equating the coefficients before the same powers of 

X - x-(x^_^i-x^)/2. This results in the following formulae for the functions yi\x), 

y^'\x)mdy^^\x) 

(x) = x-^gJ'^x3+... , (6b) 

4\x) = \li'b^ + '-Ri'b\.. . (6c) 

Here we have assumed the power expansions of the functions G(x) and R(x) in the 
form 

G{x) = G^'^ + Gf + ^ G^'^jc ' +. . . , (7a) 
R{x) = R^'^ + Rl'^x +^R^'^x^+... . (7b) 



Using such solutions is similar to employing the finite difference or finite element 
method of the corresponding order. For this reason this method with using poljmomial 
functions has no evident advantages before the widely used standard methods. 

In this paper we study another form of approximation of the solutions of the 
homogeneous equation: 

yi^ix) = Ag(^)exp(4)g(x)) , (8) 
where A^lix) and ^i2(^) are the polynomial functions of the following form: 

4:^(^) = (A))2+(Ai)2x+^(A2)«x^+... , (9a) 



4) 



'^Hx) = {ko)p + kk,)p'+.... (9b) 
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At the moment we will not pay attention to the form of the specific solution which 
should depend on the form of the function R(x) in equation (1). For example, if the 
function R(x) is sufficiently slowly varying at the segment Si , a polynomial 
approximation (7b) is proper. 

Now discuss the approximation (8), (9) in more detail. Such representation of 
the solutions makes it possible to describe rapidly oscillating or exponentially growing 
(decreasing) solutions at the segment Si even when the functions A(x) and ^>(x) vary 
slowly at this segment. (Here and below we omit the indices of the solution number 
and the number of segment). Such situation takes place in the case when Goh^>l, but 
Gih^«l. In this situation the polynomial approximation of the solutions (6a), (6b) 
leads to large errors. When Goh^«l we expect that the pol5aiomial-exponential 
solutions and the polynomial ones will behave similarly. 

First it is necessary to show that such type of solutions can fulfill the 
homogeneous equation (1) with certain accuracy. After substituting them into the 
equation we obtain 

A'\x) + 2AXx)k(x) + A(x)k\x) + [k'^(x) + G(x)\A(x) = , (10) 

where prime denotes the derivative by x. Here we have introduced the function 
k(x) = <I>'(x) . The left-hand side of equation (10) is a pol5aiomial and, therefore, to 
fulfill this equation with prescribed accuracy we have to nullify the coefficients before 



different powers of x from the zero one up to the power corresponding to the degree 
of the approximation. This results in the following equations 

(kQ +Gq + ki^AQ + IkgAi + A2 = , (11a) 
(Gi + IkQki +k2)AQ+ + <^0 + 3^1)^1 + 2^o^2 + A3 = , (1 lb) 

In practice we have to cut the series (9) at some maximum powers, jtia for A(x) and 

m$ for4>(^), sufficient to obtain the prescribed degree of approximation. To 

determine the unknown coefficients and km we require niA+m^ equations. They can 
be picked up from the set (11). Note that since the set of equations (11) is 
homogeneous, one of the coefficients Am , e.g. Ao, should be assumed to be known. 
Note that the set of equations (11) is Unear in coefficients A,„ and non- linear in 
quantities km. In general case such system of equations is difficult to solve. We can 
simplify the problem decreasing the number of equations involved from (11), which 
reduces the degree of approximation. This creates some freedom in choosing the 
quantities km . The most simple way is to determine k(x) from the equation 

k^(x) = -G(x) , (12) 
which is similar to the zero-order WKB approximation for equation (1). In this case 
we use the first niA equations from (11) to find the coefficients Am . For mA=2 and 
m$=l we can readily obtain the formulae 

^0=±V^,^l = ±^-p^, (13) 

Aq = -(aGq + 3ki)c , Ai = -IkQkiC , A2 = 3kiC , (14) 
where C is an arbitrary constant. 

Another method of obtaining the local solutions we illustrate for mA= m*=l . We 
shall use two first equations from the set (11). This leads to an additional condition for 
quantities ko and ^1: 

3ki + 4kiGo + [kQ + Gof - IkgGi = . (15) 



If we put kQ = -Gq , which is the solution of equation (15) for mA=m^=0 (zero- 
order approximation), then we obtain the quadratic equation which gives us ki (first- 
order approximation). Thus, in the framework of this scheme we have 



^0 = -^|-Go , ki=--Go 
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(16) 
(17) 



\(0 



Aq = 2k()C , Ai = -kiC . 
From (16) foUows that (^o)|'^ ""(^0)2^ but {kif'^' ^-{ki)^^' • This asymmetry of 
k{x) is in contrast with the WKB-t)^e relation (12). 

This scheme can be continued by consequent simultaneous increasing ma and m$ 
by unity. In this case equation (15) will be more complex. However, we will get a 
quadratic equation for the highest quantity A:^^ , if we substitute km values from lower 

order approximations. 

Another modification of the scheme described can be obtained using the 

requirement of symmetry of two roots of k{x), instead of the assumption k^ = -Gq . 
In this case odd and even parts of equation (15) can be separated, which yields two 
equations for ko and h . Finally, we obtain the following formulae 
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Ai 
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k^=± 
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(18) 



Aq = IkQC , Ai = -{ko + Go + ki)c 



(19) 



Further we shall call the local solution method in the form (13) and (14) as 
version 1, in the form (16) and (17) as version 2, and in the form (18) and (19) as 
version 3. 

Note that all the local solutions of the form (8) considered above tend to the 



WKB solutions in the limit 



Go' 



» 



Gl . This makes it possible to use large-step mesh 



for nimierical calculations even in the regions where WKB approximation is valid. 

When using the local solutions in the above mentioned forms one should keep in 
mind that they either diverge or degenerate at the segment where Go=0 . Besides, there 
can arise some other points where the local solutions are degenerate. This is a 



disadvantage as compared with polynomial solutions which are always suitable. The 
points of degeneration can be found from analyzing the Wronskian for the solutions. 

For version 1 of the method the condition of Wronskian nullifying at the segment is 

Go(l5Gf +64Go) = , (20) 
which yields again the condition Go=0 and, besides this, a new degeneration point that 



appears at some negative Go value. Since here 



Go' 



the point is situated in the 



region where the WKB approximation is not vaUd and the solution does not oscillate 
or vary rapidly. Thus at the segment of degeneration the polynomial solution can be 
used instead. 

For version 2, the condition of Wronskian nullifying reads Gq = . In this case 

no additional bad points arise. 

For version 3, the condition is 

Go(3Gi^+16Go) = . (21) 

Essentially, this condition is similar to the condition (19) except for numerical 
coefficients. 

3. Numerical experiments 

In this section, using numerical tests, we compare the versions of the local 
solution method proposed above with the standard finite difference method. For these 
tests we use the homogeneous Helmholtz equation (1) putting R(x)=0. In this case, 
non-trivial solutions of the equation appear owing to non-homogeneous boundary 
conditions. For the finite difference method [6] we employ the uniform mesh for which 
the finite difference scheme can be written in the following simple form 

yi+1 -'2-yi + yi-l + h^G{xi)yi = O , (22) 
where h = xi- xi_i . 

The local solution method has been realized numerically following the procedure 
described in section 2. All three versions of the considered local solutions have been 
tested. The main part of calculations was performed for the function G(x)=3x . Thus, 

the solution can be expressed as a combination of Airy functions, Ai(V3x) and 



Bi(^/3x). Besides this, the function G(x)=3x-O.06x^ has been used. The boundary 
conditions employed are the following 

As has been pointed out above, under certain conditions the local solutions in the 
exponential form can degenerate. We use two methods to overcome this difficulty. The 
first one is to use the pol5aiomial local solutions at bad segments instead. The 
degeneration can be also avoided if we shift the end points of the bad segment. In the 
calculations we use both methods. In the case of version 1 the first method has been 
used. For version 2 we have used the second method. Remind that for this version the 
only case of degeneration is Gq ~ 0. For version 3 we improve the segment with 

Gq ~ by the second method and use the first method at segments containing the 
additional points of degeneration. 

As a quantitative characteristics of the calculation accuracy we introduce the 

local error 

Sjj = yi-yex{xi) , (24) 
and the relative error 



5 = 



I 



L{yi-yex{xi)Y 



L Jex(^j) 



(25) 



where Jex(^) is the exact solution. As yex(^) we have used a solution obtained at 
much finer uniform mesh. The summation in (25) is performed over the fine mesh 
nodes. 
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Figure 2. Solution of equation (1) for G(x) = 3x, x/= -2.5, x,=1.5. 
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Figure 3. Error of numerical solution of equation (1) showed in figure 2, 
(a) - for version 3 of local solutions method, (b) - for version 2, (c) - for 
version 1 and (d) - for finite difference method. 



Figures 4a, 4b and 4c display the relative errors as a function of the number of 
mesh nodes A'^ for all three versions of the local solution method. For comparison, the 
corresponding curve for the finite difference method is shown at every plot of figure 4 
. First of all, we have to note that all the versions demonstrate the convergence to the 
exact solution in average. The average rate of convergence is proportional to 
which corresponds to the degree of approximation and is similar to what we have in 
the case of the finite difference method. At the same time, all three versions yield more 
than two orders lower error level as compared to the finite difference method on the 
same mesh. A characteristic feature of the local solution method is non-monotonous 
convergence. The oscillations of the relative error are the largest in the case of version 
1. For versions 2 and 3 these oscillations tend to vanish when A'^ increases. 

We have also tested the method proposed for a solution having much more 
oscillations in the x domain. In this case we have used another form of the G-function, 

the parabolic one: G(x) = 3x- 0.06x . For the calculations we choose x/= -2.5, 
x^45.0 . The corresponding solution is shown in figure 5. It can be reproduced by the 
finite difference method. At A^=20000 it jdelds the relative error 5= 2.110"l Version 3 
of the local solution method yields nearly the same error 5= 2.2-10"^ at A^=240. 
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Figure 4. Relative error vs. number of mesh points, (a) - for version 3, (b) 
-for version 2 and (c) - for version 1. In all three figures right curve displays 
relative error for finite difference method. 




0.0 20.0 40.0 
X 

Figure 5. Solution of equation (1) for G(x) = 3x- 0.06x , xi- -2.5, 
x,=45. 



4. Conclusions. 



In the presented paper we have analyzed and tested numerically the local 
solution method as applied to the Heknholtz equation. This method differs in principal 
from the standard methods of finite differences and finite elements. The method of 
local solutions makes it possible to use a wider class of basic fimctions approximating 
the solution, not only the polynomial ones. The most gain can be obtained when the 
basic ftinctions are close to exact solutions. In the present paper we have considered 
the basic functions of the polynomial-exponential type. The form of these functions is 
close to the WKB solutions. Owing to this we have obtained good numerical accuracy 
if the solutions are rapidly oscillating or exponentially increasing (decreasing). Besides, 
these basic functions are also able to reproduce the solutions with prescribed accuracy 
in regions of slow solution variation. This allows one to apply the method of local 
solutions to general problems. 

The numerical tests performed have demonstrated convergence of the method 
proposed. In average the convergence corresponds to the degree of approximation, 
although it is not monotonous general. The method proposed has shown the evident 
advantage before the standard finite difference method of the same order for modelling 
solutions of the Airy equation and the parabolic-cylinder-type equation. Depending on 
the solution character, the local solution method requires 10-100 times rarer mesh than 
finite difference method to obtain the same accuracy. 

References 

1. Jaeger, E.F., Batchelor, D.B., Carter, M.D., Weitzner, H.: Nucl. Fusion 30 (1990) 
505 

2. Gruber R. and Rappaz F. "Finite-Element Methods in Linear Ideal MHD" (Springer- 
Verlag, BerUn, 1985) 

3. Moiseenko V.E., Lyssoivan A.I., Pljoisnin V.V. et al, in Proc. 1996 Int.Conf. on 
Plasma Physics (ICPP 1996), Nagoya, Japan, 1996, Ed. H.Sugai and T.Hayashi, NIFS, 
Vol.2, P. 1346. 

4. BirdsaU C.K., Langdon A.B "Plasma Physics via Computer Simulation" (McGraw- 
Hill Book Company, 1985) 



5. Handbook Of Plasma Physics. Vol.1. Basic Plasma Physics, Ed. by Galeev A. A. 
and Sudan, R.N. (North-Holland Publ. Company, 1983) 

6. J.M.Ortega, W.G.Poole , Jr. "An Introduction to Numerical Methods for 
Differential Equations" (Pitman Publishing Inc., 1981) 



